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Abstract. 

We propose a simple model based on the Gnedenko limit theorem for sim- 
ulation and studies of the ordinary Levy motion, that is, a random process, 
whose increments are independent and distributed with a stable probability 
law. We use the generalized structure function for characterizing anoma- 
lous diffusion rate and propose to explore the modified Hurst method for 
empirical rescaled range analysis. We also find that the structure function 
being estimated from the ordinary Levy motion sample paths as well as the 
(ordinary) Hurst method lead to spurious "pseudo-Gaussian" relations. 

PACS number(s): 02.50.-r, 05.40. +j 

By Levy motions, or Levy processes, one designates a class of random 
functions, which are a natural generalization of the Brownian motions, and 
whose increments are stably distributed in the sense of P. Levy[|l|. Two 
important subclasses are (i) ordinary Levy motions (oLm's), which gener- 
alize the ordinary Brownian motion, or the Wiener process ||, and whose 
increments are independent, and (ii) fractional Levy motions (fLm's), which 
generalize the fractional Brownian motions (fBm's) || and have an infinite 
span of interdependence. 

The Levy random processes play an important role in different areas of 
applications, e.g., in economy [Q, biology and physiology ||, fractal and 
multifractal analysis M, problems of anomalous diffusion M etc. In this 
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paper we provide a simple method for numerical simulation of oLm's and 
discuss their scaling properties. 

At first we show the way to generate random sequence of independent 
identically distributed (i.i.d.) random variables possessing stable probability 
law. These variables play the role of increments of the oLm. 

We restrict ourselves by symmetric stable law with the stable probability 
density p a ^(x) and the characteristic function 

PaAk) =< e ikx >= exp {-D \k\ a ) (1) 

Here a is the Levy index, < a < 2, and D is a positive parameter. At 
a = 1 and 2 one has the Cauchy and the Gaussian probability laws, respec- 
tively. In other cases the symmetric stable laws are not expressed in terms 
of elementary functions. At < a < 2 they have power law asymptotic tails 

, , + a) sin(7ra/2) , . 

p a;D (x) oc D— \ 1 \ x -> ±oo (2) 
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Among the methods of generating random sequence with the given proba- 
bility law F(x) the method of inversion seems most simple and effective 
However, it is well-known fact that its effectiveness is limited by the laws 
possessing analytic expressions for F~ l , hence, the direct application of the 
method of inversion to the stable law is not expedient. In this connection, 
we exploit an important property of stable distributions. Namely, such dis- 
tributions are limiting for those of properly normalized sums of i.i.d. random 



variables [10]. To be more concrete, we generate the needed random sequence 
in two steps. At the first one we generate an "auxiliary" sequence of i.i.d. 
random variables whose distribution density F'(x) possesses asymp- 

totics having the same power law dependence as the stable density with the 
Levy index a has, see Eq.(2). However, contrary to the stable law, the func- 
tion F{x) is chosen as simple as possible in order to get analytic form of F -1 . 
For example, 

f ^(l + ixDp 1 x<0 
HX) - 1 l-[2(l + s«)]-\ x>0 (3) 

At the second step the normalized sum 
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where 

1/q 
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is estimated. According to the Gnedenko theorem on the normal attraction 
basin of the stable law lfL0|| , the distribution of the sum (4) is then converges to 



the stable law with the characteristic function (1) and D =1. It is reasonable 
to generate random variables having stable distribution with the unit D, with 
a consequent rescaling, if necessary. Repeating N times the above procedure, 
we get a sequence of i.i.d. random variables {X n {m)} , n = 1, ...,N. In the 
top of Fig.l the probability densities p(x) for the members of the sequence 
{X n (m)} (m = 30) are depicted by black points for (a) a = 1.0, and (b) 
a = 1.5. The functions p a ,i{x) obtained with the inverse Fourier transform, 
see Eq.(l), are shown by solid lines. In the bottom of Fig.l the black points 
depict asymptotics of the same probability densities in log-log scale. The solid 
lines show the asymptotics given by Eq.(2). It is seen that the Levy index can 
be estimated with the use of X n 's, which lie outside the peak located around 
x = 0. The examples presented demonstrate a good agreement between the 
probability densities for the sequences {X n } obtained with the use of the 
numerical algorithm proposed and the densities of the stable laws. 

We would like to stress that a certain merit of the proposed model is its 
simplicity. It is entirely based on classical formulation of one of the limit 
theorems and can be easily generalized for the case of asymmetric stable 
distributions. It is also allows one, after some modifications, to speed up 
the convergence to the stable law. These problems, however, ought to be 
the subject of a separate paper. We note, that two schemes were proposed 
recently, which use the combinations of random number generators |TT| and 



the family of chaotic dynamical systems with broad probability distributions 
EJ, respectively. The former method allows one to generate the sequences 



with the symmetric laws, whereas the latter allows one to generate also asym- 



metric ones. The comparison between our scheme and those of Refs. f|TT| , |12 
is beyond the scope of our paper. In any case, the proposed method can 
serve for a further constructing of non-stationary processes and studying of 
their properties. We proceed to this task below. 
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With the help of the sequence obtained, the oLm is defined by 

L a (t) = J2Xn (6) 

n=l 

(below we denote time argument as t, r for the continuous and for the discrete 
time scales as well; t, r take positive integer values in the latter case). 

In Fig. 2 the stationary sequences of independent random variables ob- 
tained with the numerical algorithm proposed are depicted by thin lines at 4 
different Levy indexes. The thick lines depict the sample paths, or the tra- 
jectories, of the oLm's. It is clearly seen that with the Levy index decreasing, 
the amplitude of the increments increases. The large sparse increments lead 
to large "jumps" (often named as "Levy flights") on the trajectory. 

Let us proceed with the properties of self-similarity of the oLm. The 
characteristic function of the oLm increments is 

(exp [ik(L a (t + t) - L a (t)]) = exp(- \k\ a r) (7) 

(D = 1 here and below). The increments of the oLm are stationary in a 
narrow sense, 

L a (ii + r)-L a (t 2 + r)=L a (ii)-L a (* 2 ) , (8) 
and self-similar with parameter 1/a, that is, for an arbitrary h > 

L a (t + r) - L a (t) £ {h- 1 ^ [L a (t + hr) - L a (t)}} , (9) 

where = implies that the two random functions have the same distribution 
functions. 

We consider two corollaries of Eqs.(7) - (9). 

1. A "1/a law" for the generalized structure function (GSF) of the oLm 
can be stated as follows: for all < \i < a the l//i-th order root of the GSF 
is defined by 

S^ir, a) = (\L a (t + r) - L a {t)\») l '» = r^V^; a), (10) 
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where 



V{fi;a) = lJ dx 2 \x 2 \' 1 J exp(-^ix 2 - I (11) 

k— oo oo J 

For the ordinary Brownian motion a = 2, and 1/2 law is the indicator 
of classical (normal) diffusion. Since r-dependence is not changed with /x 
varying, then the quantity Sl/^(r;a) at any /x less than a can serve as a 
measure of anomalous diffusion rate. We remind that the (ordinary) structure 
function is infinite for a < 2. 

We study numerically the dependence of the index s in the relation 

S^(r;a)^r s (12) 

vs fi,a. In Fig. 3 s vs a is depicted by crosses at fixed /x = 1/2. The 1/a 
curve is shown by primes. One can be convinced himself that the 1/a law 
for the GSF is well confirmed at /x smaller than the smallest Levy index in 
numerical simulation. For the comparison s vs a is depicted by black points 
for the structure function, that is, for /x = 2. It is shown that the structure 
function, being estimated from a finite sample path, lead to the spurious 
"pseudo- Gaussian" value s — 1/2. At the inset s vs /x is depicted for the 
oLm with a = 1. It is shown that s = 1 at /x < 1 , whereas with /x increasing 
the deviation from 1/a law increases. 

The "pseudo-Gaussian" r— dependence of the structure function can be 
explained by the finiteness of sample length taken into account. Indeed, let 
J max be the mode of maximum value (that is, the most probable maximum 
value) for the sequence {X n } , which consists from t terms having stable 
distribution with the Levy index a. It can be easily shown that X max oc t 1 ^, 
and for the variance we get (X 2 ) oc t 2 /" -1 . Therefore, for r smaller than 
t (which is the natural condition when estimating the structure function in 
numerical simulation or at data processing) one gets in the discrete time scale 




(L a (t + T)-L a (t)f) = W/X^ocrt 2 /"- 1 . (13) 



Thus, the square root of the structure function behaves as r 1//2 for all a's, as 
it is indeed demonstrated in Fig. 3. Furthermore, the value of the structure 



5 



function grows with the number t of the terms in the sample path growth. 
On the contrary, the GSD of the /i-th order, being estimated from a finite 
sample path at // < a, does not grow with t growth. 

2. A "1/at law" for the span of L a (t) can be stated as follows: 

R(t) = sup L a (t + r) - inf L a (t) = r 1/a R(l) (14) 

0<t<T 0<i<r 

For the ordinary Brownian motion r _1//2 i?(r) has a distribution independent 
of t|I^||. In the empirical rescaled range analysis, that is, at experimental 



data processing or in numerical simulation the span of the random process 
is divided by the standard deviation (that is, the square root of the second 
moment) for the sequence of increments, which "smoothes" the variations 
of the span on the different segments of time series WM. Such a procedure, 



called the Hurst method or the method of normalized span, is not satisfactory 
for the Levy motion because of the infinity of the theoretical value of the 
standard deviation. Therefore, we propose to modify the Hurst method by 
exploiting the 1/a -th root of the a-th moment instead of standard deviation, 
that is, 

<y a =[-Y.\ x n\ a ] (15) 

Since it has only weak logarithmic divergence with the number of terms 
in the sum increasing, then one has 




(16) 

where H = 1/a is the Hurst index for the oLm with the Levy index a, and 
the bar denotes averaging over the number of segments (having the length 
r) of the sample path. 

Fig. 4 demonstrates the application of the modified Hurst method to the 
sample path of the oLm with the Levy index a = 1. In Fig. 4a the fluctuations 
of span (thin curve) and those of <j a (thick curve) are shown for the case, 
when the total length of the sample is divided into 64 segments, each of 
t = 16 lengthwise. Below the variations of the ratio R/a a are depicted. It is 
shown that fluctuations of the ratio is much smaller than those of the span. 
This circumstance justify the use of the ratio in the empirical analysis. In 
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Fig. 4b the rescaled span vs time interval r is depicted in log-log scale by 
black points. The slope of the solid line is equal to H = 0.9. 

In Fig. 5 the Hvs a is depicted by crosses, whereas the curve 1/a is indi- 
cated by primes. The Hurst index obtained from a "traditional" ratio -R/02 
is shown by black points. It follows from the figure that a a only "smoothes" 
the variations of span, thus leading to the correct value H = 1/a. On the 
contrary, the standard deviation, being used in empirical analysis of the oLm, 
"suppresses" the variations of the span, thus giving rise to the spurious value 
H « 0.5. As in case of the structure function, this "pseudo-Gaussianity" is 
explained by the finiteness of the sample of the oLm. Indeed, for the segment 
of the length r R(r) oc r 1 ^, whereas a 2 oc r 1 /"^ 1 / 2 , and thus, -R/02 oc r 1 / 2 
for all a's. This circumstance allows us to suggest that at estimating the 
(ordinary) structure function and normalized span from experimental data 
the "Levy nature" of them can be easily masked. This, in turn, poses an 
interesting task of developing statistical methods for extracting reliable char- 
acteristics from experimental data, for which Levy statistics can be expected 
from e.g., some physical reasons. 

At the end we note that introduction of correlations into the sequences of 
i.i.d. stably distributed random variables allows us to extend the presented 
methods on the studies of fractional Levy motions. 

This work is done within the framework of the Project "Chaos-2" of 
the National Academy of Sciences of Ukraine and the Project INTAS 93- 
1194. The information support within the Project INTAS LA-96-09 is also 
acknowledged. 
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FIGURE CAPTIONS 

1. Probability densities (above) and their asymptotics (below) are indi- 
cated for the sequences of random variables generated with the use of the pro- 
posed numerical algorithm at the Levy indexes (a) a — 1 , and (b) a = 1.5. 
The probability densities and the asymptotics of the stable laws are indicated 
by solid lines. 

2. Stationary sequences (thin lines) and ordinary Levy motion trajectories 
(thick lines) at the different Levy indexes. 

3. Plots of the exponent s in Eq.(12) versus the Levy index a at ji — 1/2 
(crosses) and /i = 2 (black points). The 1/at curve is depicted by dashed line. 
At the inset s vs /i at a = 1 is shown. 

4. (a) The variations of span (thin curve), of the GSD of the cc-th order 
(thick curve) and of their ratio (below) at the different time intervals for the 
oLm with a — 1. (b) Rescaled span vs time interval in log- log scale (black 
points). Solid line has a slope H = 0.9. 

5. Plots of the Hurst exponent H vs a estimated with the use of Eq.(16) 
(crosses) and with the use of the "traditional" Hurst method (black points). 
The 1/a curve is depicted by dashed line. 
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